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ABSTRACT 

Nonlinear approximation methods such as the Zeldovich 
approximation, and more recently the frozen flow and linear potential 
approximations, are sometimes used to simulate nonlinear gravitational 
instability in the expanding Universe. We investigate the relative 
accuracy of these approximations by comparing them with the exact 
solution using second order perturbation theory We evaluate the density 
and velocity fields in these approximations to second order, and also 
determine the skewness parameter £3 = (<5 3 ) / ((S^) 2 ) 2 for each of the 
approximations again in second order. We find that S3 = 4, 3, 3.4 for 
the Zeldovich approximation, the frozen flow and the linear potential 
approximations respectively as compared to S3 = 34/7 for the exact 
solution. Our results show that, of all the approximations considered, the 
Zeldovich approximation is the most accurate in describing the weakly 
nonlinear effects of gravity. Moreover, the Zeldovich approximation is 
much closer to the exact results for matter and velocity distributions 
than the other approximations if the slope of the power spectrum of 
density perturbations is — 3 < n < —1. 

Subject headings: Cosmology, gravitational clustering, nonlinear 
approximations, non-Gaussian statistics. 
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1. Introduction 



It is now commonly agreed that all gravitationally bound objects in the Universe 
as well as its large-scale structure originated from the growth of initially small 
inhomogeneities in the expanding Universe. Generally speaking the investigation 
of non-linear evolution of non-interacting particles constituting the dark matter in 
the Universe requires long N-body computer simulations and as a result has been 
performed only for a selected subclass of initial conditions. That is one of the reasons 
why different approximation schemes have been proposed which greatly reduce the 
required numerical work and make it possible to investigate non-linear evolution for 
a much wider class of initial conditions and for longer periods of time. 

Amongst nonlinear approximation methods the following three are the most 
natural and straightforward: 

1) the Zeldovich approximation (Zeldovich 1970) and its further generalization to 
the period after caustics formation - the adhesion model (Gurbatov, Saichev & 
Shandarin 1985); 

2) the frozen flow approximation (Matarrese et al. 1992); 

3) the linear potential approximation (Brainerd et al. 1992; Bagla & Padmanabhan 
1993). 

All these approximations are really approximations in the sense that they are neither 
exact nor asymptotic to the exact solution beyond linear order (apart from some 
special degenerate cases). The difference between each of these approximations and 
the exact solution arises already in the second order of perturbation expansion. 
Thus, a natural way to see the difference between these approximations and to 
estimate their relative accuracy is to calculate their departure from the exact solution 
in this order. This is just the aim of the present paper. 
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2. 



Second Order Perturbative Calculations 



Let us consider gravitational instability in the spatially flat matter dominated 
FRW Universe before the formation of caustics. In such a Universe the scale factor 
grows as a(t) oc t 2 ^ 3 and the background matter density decreases according to 
p = 1/Q-rcGt 2 . The equations describing this process have the following form in the 
Newtonian approximation ( see, e.g. Peebles 1980): 



where 5 = (p - p )/p o . 

The velocity u is irrotational, therefore it is possible to introduce a velocity 
potential V so that u = —W/a. Then Eq. (3) may be substituted by its first 
integral: 



A$ = 4irGa 2 p 5; 

5 + -div((l + 5)u) = 0; 
a 

(au)' + («V)u= -V$ 



(1) 



(2) 



(3) 



U = $ + 4(W) 2 , 



(4) 



and the first two equations acquire the form: 




(5) 



a 2 5 = (1 + 5)AV + V5W. 



(6) 



Let us expand all quantities into series in powers of an initial density 
enhancement: 5 = 5^ + 5^ + and the same for $ and V. As a function of 
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time, this is an expansion in powers of t 2 /a 2 oc t 2//3 (since 5^ oc t~^). The first 
order solution (the linear approximation) is the same for all three abovementioned 
approximation schemes: 



It 2 

$ (1) = 0o; v^ = U; 5(1) = |^>o (7) 

where (f>o(r) is the initial gravitational potential. The value and statistical properties 
of 4>o (r) are completely arbitrary for classical cosmology (apart from the trivial 
condition of smallness: |0o(^*)| <C 1)- On the other hand, any theory of the Early 
Universe should produce some predictions for the properties of fa. In particular, the 
inflationary scenario of the Early Universe predicts that O is a gaussian stochastic 
quantity with zero average and dispersion for the Fourier transform (f>o(k) having an 
approximately k~ 3 dependence (in the simplest versions of this scenario and before 
the multiplication of (f>o(k) by a transfer function). 

The form of the linear approximation (7) naturally leads to three approximation 
schemes which arise as a result of imposing by hand some of the relations valid in 
this approximation on a fully non-linear solution. Of course, having introduced one 
new relation, one has to abandon one of the previous equations. That abandoned 
equation is chosen to be the Poisson equation (1) or (5). Thus, in all these 
approximations we neglect the self-gravity of inhomogeneities. As a result, we have 
the following relations for the full series: 

1) V — &t - the Zeldovich approximation (further denoted by ZA and by the subscript 

zy, 

2) V = (pot - the frozen flow approximation (FF, subscript /); 

3) $ = 0o - the linear potential approximation (LP, subscript p). 

Our way of introducing the Zeldovich approximation is different from that commonly 
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found in the literature. However, it is straightforward to check that the formula 
V = provides both the necessary and sufficient conditions for conventional 
expressions of the Zeldovich approximation, since, (a) it can directly be derived 
from those latter expressions, and (b) by inserting V = $t in Eq. (4) and then 
solving Eq. (4) using a trivial change of variables (V = a 2 aV), we arrive at the 
Zeldovich - Bernoullie equation commonly used to describe the evolution of the 
velocity potential in the Zeldovich approximations (see eg. Kofman 1991) which is 
given as follows 

oa 2 

Let us now consider the second order terms. 



1) The Exact second order solution 



A*' 2 ' = |^ 2 >; (8) 
a 2 5^ = + S^AVM + W (1) W (1) ; (9) 

1/(2) = $(2) + J_ (W«) 2 . (10) 
2a 2 

After excluding $( 2 ) and V < - 2 ' ) from these equations, we get an equation for 5^: 



5 (2) + p (2) - J^ {2) = IJ^wavw + VSMVVM) + ((V\/«) 2 ) . (11) 

Note that the left hand side of this equation is the same as in the first order. The 
solution of (8-11) is (we consider the growing mode only): 
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3f : 



P(f) = (A0 O ) 2 + V0oV(A0 o ) = V(V0oA0 o ), 



Q(f) = (V0 O ) 2 . 



(12) 



Of course, this solution is well known (see, e.g. Peebles 1980). The average value of 
is zero because it has the form of divergence (the same is true in all orders). The 
expression for 5^ is local while the appearance of the inverse Laplacian A -1 in the 
expressions for <3>( 2 ) and shows that they are non-local, as is the velocity 
(but divu^ is local). Non-local terms in the expansion of 5 begin from 5^, thus, 
the first nonlinear corrections to the power spectrum and the density correlation 
function are non-local, too. 

2) The Zeldovich approximation 



4 2) = ^; ^i 2) = 4 2) + ^(v^)) 2 , (is) 

the third equation is the same as Eq.(9). The solution is: 



^ = i!i (2P+aq); v ^ = ^ = h Q - (14) 

Note that in the case of one-dimensional plane-symmetric motion O = <Po(x), 
P = §AQ and the above terms coincide with the second-order terms for the exact 
solution (12). This is a consequence of the fact that the Zeldovich approximation 
is actually an exact solution of Eqs. (4 - 6) in the case of one-dimensional motion 
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before caustic formation (Shandarin & Zeldovich 1989). 



3) The frozen flow approximation 



= *. *? = -h {vvWf = ~& Q - (15) 

the third equation being the same as Eq.(9). The only quantity that remains to be 

(2) 

found is 5f . From the above equations it follows that: 



5f = *Lp. (16) 
f 8a 4 K ' 



4) The linear potential approximation 



$f = 0; V?> = ^(W«) 2 , (17) 

the third equation still being the same as Eq.(9). Solutions for the remaining 
quantities are 

W-mGW-m?^*"*- (18) 

Note that 5^ = 0.75^ 2 \ Thus, the second order correction in the linear potential 
approximation has the same spatial structure as the exact solution but its value is 
30% smaller than that of the exact solution. 
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3. Comparison of Approximations 

3.1. Density Perturbations 

Comparing second order terms in the density perturbation in the different 
approximation schemes, we see that the P terms are the same in all of them (and 
different from the exact solution). The difference between 5^,5^ and 5^ arises 
essentially because of the different numerical coefficients in front of Q in each of 
these approximations. In addition we would like to point out that there are no 
nonlocal terms in 5^ in any of these approximations. 

Let us now consider the difference A between the approximate solutions and 
the exact one in the second order, and also calculate the expected variances of A 
assuming that the initial potential (poif) is a gaussian stochastic quantity with zero 
average and an isotropic power spectrum. We have 



A P = 6® - <5< 2 > = _^( P + I A Q) = -0.3^. (19) 

We introduce the notations: a\ = ((V</> ) 2 ); a\ = ((A0 O ) 2 ); °l = ((V(A0 o )) 2 ). 
In the linear approximation, a\ is proportional to the velocity dispersion 
(a 2 = ((u^) 2 ) = ^-o" 2 ), of - to the dispersion of density perturbations 
(a 2 = ((5^) 2 ) = ^cr 2 ). of < a±a 3 with the equality being achieved in the case when 
the Fourier spectrum is proportional to S(k — k ) only. Using the useful relations 
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15 z 3 

after lengthy but straightforward calculations we get 

(A|) = 



I C>2\ ^ 4 | 1 2 2 

(-P ) = 3^2 + g^i^; 
(PAQ) = 2cr 2 4 + j^V 2 ; 
((Ag) 2 ) = ^a 2 4 + ^ 2 4 (20) 



2 /851 4 49 2 2 \ 2 



where A = ||^. To obtain a relative accuracy with respect to the second order 
term in the exact solution, these results should be divided by (d^ 2 ) = ^(A 2 ). 
The form of the fractional error F$ = (A 2 )/(<5 < - 2 - )2 ) for different approximations can 
be expressed as Fs = (cry 2 + 5)/(c7 2 + 1) where 7 = <7 2 /<7i<73. The values of a, b, 
c for the three approximations considered by us are (0.15, 0.0, 6.25) for ZA; (1.1, 
0.25, 6.25) for FF; and (0.0, 0.09, 0.0) for LP respectively. From the plot (Fig.l) 
of Fs as a function of 7, it is clear that the Zeldovich approximation is always 
better than either FF or LP. In particular, for a 5-like power spectrum (7 = 1) the 
fractional error F s is 0.0246, for ZA, 0.187, for FF and 0.090 for LP. In the reverse 
case of a very extended spectrum 7 <C 1 the fractional errors have the asymptotic 
forms 0.157 2 ; 0.25; 0.09 for ZA, FF, LP respectively. (For a power-law spectrum 
(— )| oc k 11 , 7 <C 1 if the spectral index lies in the range — 5 < n < — 1). Note that 
in the latter case the Zeldovich approximation is much closer to the exact solution 
than the other two approximations since ((A|) <C (A 2 J). 

These results may also be used to compare the value of the skewness parameter 
S3 = (5 3 )/o"p which arises in each of the above approximations with that obtained 
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in the exact solution (Peebles 1980, Grinstein & Wise 1987, Bouchet et al. 1992). 
To first order in perturbation theory (5 3 ) = ((<^) 3 ) = 0. In second order, however, 
(S 3 ) = 3((5^) 2 5^), so that S 3 = f « 4.86; 4; 3; 3.4 for the exact solution and 
for the Zeldovich, frozen-flow and linear potential approximations respectively. 
All three approximations produce a low value for the skewness but the Zeldovich 
approximation is the closest to the correct answer once more ( being accurate to 
within 20% ). It is also interesting that in all four cases the skewness does not 
depend upon the form of the initial spectrum. 



3.2. Peculiar Velocities 

We now proceed to calculate the error in the peculiar velocity for each of the 
approximations considered earlier. Let = —u^ denote the difference between 
the second order velocity field in the ZA and in the exact second order analysis 
(similarly for FF and LP). Let us also define M 2 = ((V(A- 1 P)) 2 ) = -(PA^P) 
(this will be the only non-local term in the answer). Using the useful relations 



(V(A^P)Vg) = -(PQ) = \a\a\- 



(VgVQ) = ^a\al- (22) 



we get in the second order: 



<^ 2 > = B\M 2 - \a\al)- 
(Uf 2 } = B 2 (M 2 + ^ala 2 2 ); 

(U^) = B 2 (M 2 + ^-ay 2 ), (23) 
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where B = In order to obtain a relative accuracy, each of these expressions 

should be divided by (u^ 2 ). Since uf^ = 0, Jjf^ = —u^ 2 \ therefore this is equivalent 

—* (2)2 

to dividing by (Uj ' ). From Eq. (23) it is clear that of the three approximations the 
Zeldovich approximation is always closest to the exact solution. 

The detailed expression for M 2 in terms of Fourier components of the initial 
gravitational potential ((f>o(r) = (2n)~ 3/2 J d z k$^e ik? \ = 4> 2 (k)5^(k - k!) 

where k — \k\ and 5^ is now the 3D delta- function) is quite complicated (see 
appendix): 



I roo roc 

M ' = 32Wo kl<l)2 ^ dkl I fc 20 2 (^) 

( 2kMk\ + 4kjk 2 2 - k\) + {k\ - k\f In ,f 1 + M dk 2 (24) 
V " " \ki-k 2 \) 

but if the physical wavelengh 2na/k making the main contribution to o\ is much 
larger than corresponding lengths for a 2 and <t 3 , then the integrals in Eq. (24) 
decouple and M 2 ps \^\<t 2 . In this case, the Zeldovich approximation is much closer 
to the exact solution than the other two. For a power-law spectrum, this happens if 
— 3 < n < —1 (certainly, a cut-off at both some large and small scales is implicitly 
assumed). On the other hand, M 2 = afa 2 in the opposite case of a S-\ike isotropic 
power spectrum, so that in this case the relative accuracy of the approximations is 
0.27; 1.0; 0.53 for ZA, FF and LP respectively. 

Juszkiewicz et al. have recently suggested that moments of the dimensionless 
velocity divergence 9 = a~ l divu = -(3t/2a 2 )AV (chosen so that (# (1)2 ) = (5 (1)2 )) 
may be useful statistical quantities of study in the weakly nonlinear regime 
(Juszkiewicz et al. 1993). The value of the skewness parameter T 3 = (# 3 )/((# 2 )) 2 for 
6 can be determined from the results obtained in the previous section. We find that 
in second order T 3 = 3((0^) 2 6^) / ((6^) 2 ) 2 = (-f w -3.71, - 2, 0, - 0.8) for the 
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exact solution and for the Zeldovich, frozen flow and linear potential approximations 
respectively. We find that once more the ZA is the most accurate of the three 
approximations, although the accuracy of all approximations worsens in this case. 

One can also calculate the fractional error F e for just as we had done for 5 
earlier. We define the difference D z = 0^ — 0^ (similarly for Df and D p ) and 
construct fractional quantities F e = (D 2 ) / (0 ( - 2,S)2 ) (i = Z,f,p). As in the case of Fg, 
F e also has the general form F e = (a^y 2 + 6) /(cry 2 + 1). Where the coefficients a,b,c 
have the values (0.59, 0.0, 3.47) for ZA, (0.0, 1.0, 0.0) for FF, and (2.83, 0.75, 3.47) 
for LP respectively. Fg is shown as a function of 7 for the different approximations 
in figure 2. It is clear that as found earlier in the case of density perturbations, the 
Zeldovich approximation performs better than the other two approximations for all 
values of 7. 

It is also possible to calculate cross correlations between the two fields 5 
and 0. Using results derived earlier one can easily show that in second order 
X 12 = (50 2 )/a* = l(S 3 - 2T 3 ) and X 21 = {5 2 0)/a 4 p = |(T 3 - 2S 3 ). 

4. Discussion 

Our results show that the Zeldovich approximation which is the simplest 
and computationally the most cost-effective of the three approximation methods 
considered by us in this paper, is also more accurate on an average than either the 
frozen flow or the linear potential approximation when studied to second order in 
perturbation theory. This might suggest that the Zeldovich approximation is a better 
tool than either of the other two approximations with which to study overdensities 
in the weakly nonlinear regime and also to probe the dynamics of underdense regions 
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such as voids. Efforts to compare the different nonlinear approximation methods 
discussed in the present paper with N-body simulations in the strongly nonlinear 
regime are presently in progress (Sathyaprakash et al. 1993). 

In addition, the Zeldovich approximation appears to be much closer to exact 
density and velocity distributions than the other approximations if the slope of 
the density power spectrum is — 3 < n < — 1 (just the range for which we expect 
the prolonged existence of a network structure qualitatively described by the 
adhesion model after caustic formation). This also means that the deviation of 
the Zeldovich solution from the exact one in the second order is much less on an 
average than the second-order term itself. This result may be generalized to higher 
orders of perturbation theory in Eulerian space as well, i.e., (A^ 2 ) <C (5^ 2 ) and 
(U^ 2 ) <C (u^ 2 ). However, this does not mean that the Zeldovich approximation 
exactly describes a non-linear evolution even in this case because this closeness 
originates from the fact that the effect of point displacement (accurately taken into 
account by the Zeldovich approximation) is more important for averaged values than 
effects of non-linearity for this range of slopes. 

Indeed, leading terms in all orders of perturbation theory having the largest 
power of a large-scale velocity may be summed with the result (the same one for 
both the exact solution and the Zeldovich approximation, see also Shandarin 1993): 



5 = 5^(q,t); 
u = u {1 \q,t); 




(V0 o (r)) 2 ; 



(25) 
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This quasi-linear solution which is actually the linear solution in Lagrangian space 
correctly accounts for a shift of phases produced by the long-wave part of the 
perturbation spectrum but does not describe genuinly non-linear effects which are 
mainly due to the short-wave part (the substitution of f by q in the last terms of the 
expressions for <3> and q that would exactly reproduce the Zeldovich approximation 
for these quantities exceeds the accuracy with which Eq. (25) is derived). That is 
why, for instance, the error of the Zeldovich approximation in determining S3 is not 
too small. 

In a companion paper (Munshi et.al 1993) we obtained values of higher moments 
of the distributions (i.e. S 4 ,S 5 , .. and T 4 ,T 5 , ..) as well. Then this are used to 
obtain P(S) vs 5, 5 vs 9 relations and void probability distribution function and 
related quantities. 
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APPENDIX 
DERIVATION OF EXPRESSION FOR M 2 



-16- 



We derive the formula for M 2 used in section 3.2. By definition M 2 = — (PA~ L P) 
It is useful to perform the analysis in k space. Decomposing the gravitational 
potential in k space we get 

Mrl = J <P(k)e^d 3 k; (Al) 

P = V(V0 o A0o) = "V^p / d 3 h J d^ikkUik)^)^^^ 

' r d 3 h [ d 3 k 2 h(h + h)kl<P(h)^k 2 )e l{A+ ^] (A2) 



(27T) 

A-'P = -J- [ d 3 k, [ d 3 k 2 ^i+^M ^O^e^^; (A3) 

(A) J J (h + k 2 ) 2 

V(A- 1 P) = / d% [ d% {kl + Wfo + k ^ ^mh)e^ + ^-, (AA) 

(A) J J (ki + k 2 ) 2 

M 2 = J d 3 hj d 3 k 2 J d 3 k 3 J rf3 A;4e i(fc 1 +fe 2 )r- l (fe 3 +fc4)r 

(k 2 + hh)(kl + hh)k^kU(h + h)(h + h)) 

The only non-zero average values are: 
(0(£i)0*(£ 2 )> = 5 3 (h - hWih); (0(^)0(4)) = 5 3 (h + k 2 )<p 2 (h), (A6) 

— * — * 

the last expression follows from the reality condition (j>(—k) = 4>*{k). So, 
(kl + hhfkj ^ (kl + hk^kl + hk^kjk 2 ] , 3 , 

— =; =; 1 =; =; U K 2 

(h + h) 2 (h + k 2 ) 2 J ' 

= ToW ^ I d ' k2 A W 2 (**) {k )r ^rt 1 ^ + ^(fc? + k 2 2 )) 

( 27r ) J J (fci + k 2 ) 2 

47T-27T [°° hi , 2(h s,, [°° U 5,2, U ^ U f 1 A (h + hz^hh + Z{k 2 + k 2 )) 

Wj / w j q k 2 Hh)dk 2 jj z > 

(A7) 
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where z is the cosine of the angle between k\ and fc 2 . Performing the z integration 
using standard textbook formulas one gets the expression for M 2 

; /*oo roc 

(lk x k 2 {k\ + 4A; 2 A; 2 - A; 2 4 ) + {k\ - A; 2 ) 3 In .f 1 + dk 2 (A8) 
which we have used in the text. 
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Figure Captions 

Fig. 1 : The fractional error F$ = (A 2 )/ (5^ 2 ) in 5 for different approximations 
is shown as a function of 7 = of /<7ia" 3 . The solid, dashed, and dotted lines correspond 
to the Zeldovich, frozen flow and linear potential approximations respectively. 

Fig. 2 : The Fractional error Fq = (D 2 ) / (9^ 2 ) in 9 for different approximations 
is shown as a function of 7 = o\lo\o?,. The solid, dashed, and dotted lines correspond 
to the Zeldovich, frozen flow and linear potential approximations respectively. 



